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ABSTRACT 

We present results of nonlinear numerical simulations of gravity wave driven shear flow 
oscillations in the equatorial plane of the solar radiative interior. These results show 
that many of the assumptions of quasi-linear theory are not valid. When only two 
waves are forced (prograde and retrograde) oscillatory mean flow is maintained; but 
critical layers often form and are dynamically important. When a spectrum of waves is 
forced, the non-linear wave- wave interactions are dynamically important, often acting 
to decrease the maintenance of a mean flow. The (in) coherence of such wave- wave 
interactions must be taken into account when describing wave driven mean flows. 



1 INTRODUCTION 



Explaining the solar internal rotation remains an outstand- 
ing problem in stellar physics. Seismic observations indicate 
that the entire convection zone of the Sun rotates differen- 
tially, with the equator spinning faster than the poles, simi- 



lar to what is observed at the solar surface ( Thompson et al 



1996 2003 1. This differential rotation persists to the base of 



the solar convection zone, below which the angular velocity 
transitions to approximately uniform rotation. The transi- 
tion between the differential rotation of the solar convection 
zone and the uniform rotation of the radiative interior oc- 
curs in an unresolved layer known as the solar tachocline 
( Spiegel fc Zahn|1992 |. 

The differential rotation of the convection zone has been 
modeled in three dimensional simulations of the solar con- 



vection zone (Miesch et al. 20061, including the effect of 



an imposed subadiabatic tachocline as proposed by Rempel| 
(20051. However, in order to understand the maintenance of 



the tachocline, one needs to understand the maintenance of 
the uniform rotation of the solar radiative interior. 

Currently, there are two main theories for the uniform 
rotation of the radiative interior: one magnetohydrodynamic 
and one purely hydrodynamic in nature. The former ar- 
gues that a primordial poloidal magnetic field confined to 
the radiative interior could enforce uniform rotation via the 



Lorentz force as described by Ferraro's isorotation law ( Fer- 
raro|1937l. This model is able to maintain uniform rotation 



only if the field lines are completely confined to the radia- 
tive interior | MacGregor & Charbonneau'1999') . If field lines 



cross the convective-radiative interface the differential rota- 
tion of the convection zone is communicated into the radia- 
tive interior, contrary to helioseismic inferences. It has been 
suggested that meridional circulation driven in the convec- 
tion zone might penetrate the tachocline and keep the field 



els have success with this approach ( Kitchatinov & Rudiger 



20061, other models fail (Brun & Zahn 20061. It appears 



that the confinement and the interior angular velocity de- 



pend sensitively on the boundary conditions imposed (Ga 



raud & Brummell 20071. It remains unclear if a magnetic 



field can effectively render uniform rotation in the radiative 
interior. 

The hydrodynamic model utilizes gravity waves to ex- 
tract angular momentum from the solar radiative interior. 
Internal gravity waves (IGW) are thought to be excited in 
the stable interior by convective downwellings hitting the 
convective-radiative interface. Such waves can cause long 
range angular momentum redistribution. Whereas earlier 
models failed to recognize the "anti-diffusive" nature of grav- 
ity waves ( Kumar fc Quataert|1997 Zahn et al.|1997 1 , more 
recent models incorporate this behaviour into the theory 
( Kumar et al.|1999[ [Talon et al.|2002[ ), which can be briefiy 
explained as follows. Internal gravity waves generated at 
the convective-radiative interface propagate into the solar 
tachocline and set up shear layer oscillations (SLO), simi- 
lar in nature to the Plumb & McEwan experiment (Plumb 



& McEwan 19781. This time-dependent, depth-dependent 



confined (Gough & Mclntyre 1998]). Whereas some mod- 



(but spherically-averaged) mean zonal flow is assumed to 
have a slightly stronger prograde sense than retrograde (on 
average), a consequence of the continual spin down of the 
Sun's convection zone by the magnetic torque arising from 
the solar wind outflow. Because of the Doppler effect, pro- 
grade waves are preferentially dissipated in the upper part 
of the radiative region, allowing predominantly retrograde 
waves (transporting negative angular momentum) to prop- 
agate into the deep radiative interior. It has been proposed 
that the solar radiative interior is thus spun down, tending 
ultimately toward a state of near uniform rotation ( [Talon 
et al.|2002| ). 

Such wave driven shear flow oscillations have been ob- 
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served in laboratory experiments ( Plumb fc McEwan|1978 l, 
in the Earth's atmosphere (Baldwin et al. 20011 and in 
Jupiter's atmosphere ( Leovy et al.|1991[ ). In addition, there 
is some evidence for oscillatory flows (the source and per- 
sistence of which is unknown) in the solar tachocline ( Howe 
et al. 20001. Therefore, it is not unreasonable to postulate 



such flows in the solar radiative interior. There are, however, 
several issues with this model. First, it is unclear that such 
oscillations could develop, at least in the upper part of the 
tachocline, due to the constant bombardment by overshoot- 
ing plumes which likely disrupts the wave-mean flow inter- 
action. Second, the development of oscillatory flow depends 
sensitively on the wave spectrum and amplitudes assumed, 
which is, at best, poorly understood. In fact, the model de- 
scribed above completely neglects the wave spectrum that 
would be produced by overshooting plumes. Furthermore, 
the model does not address the effect of a latitudinally vary- 
ing differential rotation on the wave generation, which may 
play a large role in the nature of waves which propagate 



into the deep interior (Fritts et al. 19981. Finally, even if 



oscillations could develop in the simplified sense described, 
it is unclear how Doppler filtering could produce uniform 
rotation in the deep radiative interior. Even if only retro- 
grade waves are allowed passage into the deep interior, as 
soon as one of them dissipates a radial differential rotation 
is established. Subsequent waves propagating into this dif- 
ferentially rotating region will be Doppler shifted and the 
radial differential rotation enhanced. 

The study of internal gravity waves in astrophysical set- 
tings can take advantage of the vast literature on the subject 
in the atmospheric and oceanic communities. Much of the re- 
search on angular momentum transport by IGW was spurred 
by the observation of the Quasi-Biennial Oscillation (QBO), 
which manifests itself as alternating easterly and westerly 
winds in the equatorial stratosphere. While the QBO was 
first observed in 1960 ( Reed et al.||1961 Ebdon & Veryard 
19611, the first satisfactory theoretical explanation awaited 



Lindzen & Holton (19681 and Holton & Lindzen (19721 



While different flavors of this original theory exist, the fun- 
damental theoretical explanation has persisted and has been 
borne out in experiments such as the remarkable Plumb & 
McEwan (1978) experiment. However, despite the recogni- 
tion of the physical mechanism which drives the QBO and 
significant observational details, simulating it from first prin- 
ciples has been challenging. The spectrum of waves needed 
to force the oscillation is still a topic of research ( Baldwin 



et al. 2001 Dunkerton 19971. Not all Global Climate Mod- 



els (GCM) are able to produce a realistic QBO; these sim- 
ulations depend sensitively on the convective parametriza- 



tion, resolution and diffusivities (Baldwin et al. 20011. So, 



whereas the basic properties and nature of the QBO in the 
Earth's atmosphere are understood, the details are not. 

In general, most theoretical studies of wave driven mean 
flows have been restricted to a quasi-linear formulation, a 
perturbation approach that accounts for the effects of weakly 



non-linear waves on the mean flow dynamics (Lindzen & 



HoltonI 1968) [Holton fc Lindzen|1972l|Dunkerton||1997||K~ 



mar et al.||1999| |Kim fc MacGregor||2001[ ). Clearly, this ap- 
proach limits the wave dynamics and breadth of solutions 
that can be obtained. The fully nonlinear numerical simu- 
lations described in the present paper, do not employ the 
approximations and parameterizations of the quasi-linear 



models in order to better understand wave driven mean flows 
and, in particular, those interactions under solar-like condi- 
tions. These numerical simulations are conducted in an at- 
tempt to shed light on how such waves, and the mean flows 
they produce, can contribute to the rotational profile of the 
solar radiative interior. 

The remainder of the paper is organized as follows. In 
section 2, we describe our numerical technique and in section 
3, we briefly discuss the basics of wave- mean flow interac- 
tion. In section 4, we present our results. We discuss the 
major conclusions and their implications in section 5. 



2 NUMERICAL MODEL 

2.1 Equations 

The Navier-Stokes equations are solved for a perfect gas 
within the anelastic approximation. Our model is two di- 
mensional (2D), representing only the equatorial plane. The 
anelastic approximation is appropriate when the flow speed 
is much smaller than the sound speed in the medium and 
when thermodynamic perturbations are small compared to 
the reference state; both of these constraints are well sat- 
isfied in the Sun's radiative interior. However, unlike the 
Boussinesq approximation, the anelastic approximation al- 
lows a stratification in the reference state density and tem- 
perature. The following equations are solved for the fluid 
flow relative to the rotating frame of reference and thermo- 
dynamic perturbations relative to a hydrostatic reference 
state: 



V • pi? = 



dv 



+ {v-\/)v = - VP - Cgf + 2(u X n)-|- 



dT 



+ {v- V)T = -vr{^ - (7 - l)Thp)+ 

BT 
(7 - l)ThpVr + 7K[V'r + {hp + h^)-Q^] 



(1) 



(2) 



(3) 



In these equations overbars represent reference state 
variables which are functions of radius only; the variables 
lacking overbars are the perturbations, which are functions 
of time, radius (r) and longitude {(p) . We assume no flow or 
gradients in latitude. Equation (1) is the anelastic version 
of the continuity equation. In equation (2), the momentum 
equation, v is fluid velocity (with components v^ and Vr), P 
is the reduced pressure (equal to the pressure perturbation p 
divided by the reference state density ]o), C is the codensity 
{T/T + p/ {gpT)dT/dr) |Braginsky fc Robertslig g?: 'Rogers] 
|fc Glatzniaier^^2006j ), T is the temperature perturbation, g 
is gravitational acceleration and V is the viscous diffusivity. 



ity 



the rotating frame, is a constant and in the z directio: 



^ Models with differential rotation imposed have been run, but 
for any reasonable differential rotation there is no effect on the 
wave dynamics because <5f2 is much smaller than the wave fre- 
quencies explored. 
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Equation (3) is the heat equation written in terms of 
temperature. The first term on the right hand side is the 
radial velocity times the prescribed super- or subadiabatic- 
ity, dictating the convective stability as a function of ra- 
dius. dT/dr is the reference state temperature gradient and 
(7 — l)Thp is the adiabatic temperature gradient for a per- 
fect gas, where 7 is the ratio of specific heats Cp/cv and hp 
is the inverse density scale height, dlnp/dr. The thermody- 
namic diffusivity, k, is depth dependent and h^ is dlnTc/dr. 

For numerical convenience we solve equations (1) and 
(2) using a vorticity-streamfunction formulation. Taking the 
curl of the momentum equation and imposing mass conser- 
vation we get the vorticity equation: 



In 



+ {v-V)uj = (2Vi^Uj)hpVr- 



Tr d(j) 



1 dTdp 
pTr dr d(t> 



+I7V^cj(4) 



where uj is the z component of vorticity, V x u. For simplic- 
ity we have neglected the additional viscous terms due to 
the radial gradient of the viscous diffusivity. Equation (4) 
is solved together with an equation relating the streamfunc- 
tion, -i/), defined as p« = V x i/)z, and the vorticity: 



- U)p = 



aV 



+ {--hp) 



dr 



+ 



(5) 



1 d'^ip 

This formulation automatically enforces equation (1). 
The independent variables directly calculated are then the 
temperature perturbation, T, the vorticity, u), and the 
streamfunction, 1/). 

2.2 Numerical Technique 

Equations (l)-(3) describe internal gravity wave dynamics in 
the stable solar interior, neglecting the influence of a mag- 
netic field. The equations are solved in two dimensions in 
cylindrical coordinates representing the equatorial plane of 
the Sun. The computational domain extends from O.OI-R0 to 
0.71-Rq, representing only the Sun's stable radiative interior. 
The radially dependent reference state variables (density, 
temperature, subadiabaticity and gravity) are taken from a 
21st-order polynomial fit to the one dimensional standard 



solar model Model S ( J.Christensen-Dalsgaard et al. 1991 1 



The equations are solved using a Fourier decomposition in 
the longitudinal direction and a finite difference method in 
the radial direction. The solutions are time advanced using 
a second order Adams Bashforth method for the nonlinear 
terms and an implicit Crank-Nicolson method for all the 
linear terms. A spectral transform method is employed to 
compute the nonlinear terms in grid space each time step. 
With a few exceptions, the majority of the models have a 
spatial resolution of 1000 (nonuniform) radial grid points x 
512 longitudinal grid points. The model is parallelized us- 
ing Message Passing Interface (MPI). Boundary conditions 
imposed on the bottom boundary of the domain are stress- 
free, impermeable and isothermal. The conditions on the top 
are stress free and isothermal. Waves are driven at the top 
boundary via a prescribed time-dependent permeable top 
boundary condition. 



maier|2006 l. However, in these models the spectrum of mo- 
tion in the radiative interior is complicated because of the 
turbulent convective overshoot, making a sensitivity study of 
fundamental wave dynamics virtually impossible. Therefore, 
unlike previous models, here we exclude the convection zone 
and drive the gravity waves artificially at the top boundary 
of the model, which represents to the top of the radiative re- 
gion. The simplest way to produce a wave driven mean fiow 
is to excite a prograde wave and a retrograde wave with the 
same frequency, wavenumber and amplitude. These waves 
would combine to form a standing wave in the model's ro- 
tating frame of reference if there were no mean flow relative 
to this frame. 

We employ the following procedure. The independent 
variables, vorticity, streamfunction and temperature, are ex- 
panded in a Fourier series in longitude, <j), which is in radi- 
ans. 

M 

uj{r,(j),t) = y^((^^(r,i)cos(m0) + uji"^(r,t)sm{m,(j))) (6) 



'ip(r,cl>,t) 



M 

E 



(V'm(^ t) cos{m4)) + Tp^{r, t) sm{m(t>)) (7) 



T(r, 0, t) = ^ (r;^(r, t) cos(m0) + T^^ir, t) sin(m<?!.)) (8) 

m=0 

where m is the horizontal wavenumber. Separate equations 
are solved for the sine and cosine coefficients of these vari- 
ables. 

A wave is excited at the top boundary as a radial veloc- 
ity by prescribing a longitudinal and time dependent stream- 
function there. Most models are run with amplitude con- 
stant in spectral space. Thus, equal amplitude means the 
waves have an equivalent amplitude ■)/'■ This means that the 
waves have slightly different amplitudes in grid space (since 
Vr ~ mtp/pr). Variance of amplitude in grid and spectral 
space will be discussed further in section 4. We further spec- 
trally decompose ^ in frequency space 

N M 

ipir, 0, i) = 2J 2_/ (^"^>" (*") cos(m0) cos(ncroi) 



n — Tn^O 



+'4>m,n{f') cos{m<j>) sin(ncroi) 
+'4'm,ni''') sm{m4>) cos{naot) 
+V'm,Ti(^) sin{m4>) sin(n(Toi)) 



(9) 



Here, ao is the lowest frequency (in radians/second) con- 
sidered. After some algebra (D.O. Gough, private commu- 
nication) one can show that, for a positive m, the ampli- 
tudes of a prograde wave {m4> — naot) and a retrograde wave 
inuj) + nuot) are, respectively, 



P = 2^i'4^^,n +^n 



+ {r„ 



^n 



R 



iurr, 



'rm,n} "T~ \'i-^m,n ~t~ Ym,n) 



(10) 

(11) 



2.3 Wave Driving 

More realistic models of the solar radiative interior also sim- 



ulate the dynamics of the convection zone ( Rogers & Glatz- 
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The simplest way to force a wave driven oscillation is to 
force prograde and retrograde waves with the same ampli- 
tude. By inspection of equations (lO)-(ll) one can see that 
a way to do this would be to force only one component of 
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the four i/'^,„,'(/'^,„,'(/'^,„,'(/'^,^ we choose -(/j^^^. That is, 
we force one or more wavenumber modes tp^{r,t) on the 
top boundary with one or more frequency modes cos(n(Toi) 
relative to the rotating frame freference. 



3 FUNDAMENTALS OF WAVE DRIVEN 
SHEAR FLOWS 

The dynamics of wave driven mean flows has been studied 
extensively ( |Holton|2004l |Lindzen|1990| |Plumb|1977[ ). Here 
we review some of the basics. When a single propagating 
wavaj is attenuated it transfers angular momentum to the 
mean flow ( Eliassen fc Palm|1960 1 . The forcing of the mean 
flow by waves depends on the wave attenuation, which de- 
pends on the frequency of the waves relative to the mean 
flow, which in turn depends on the shear in the mean flow. 
Therefore, wave-mean flow dynamics is highly nonlinear. 

One can appreciate how angular momentum is trans- 
ferred from waves to the mean flow by considering the 
longitudinally-averaged horizontal component of the mo- 
mentum equation for a simple constant density fluid ( |Holton| 
2004||Lindzen|1990[ ). 



dU l dr< u'w' > 
dt r dr 



8Pu_ lac/ 

Q^2 J. Qj. 



(12) 



where u' is the azimuthal velocity, v^, and w' is the radial 
velocity, Vr- This equation shows that the mean zonal flow, 
U, is driven by the divergence of the horizontally-averaged 
Reynolds stress (HARS) and is smoothed by viscous dis- 
sipation. A simple way to understand this transport is to 
consider the product of two waves, one with wavenumber 
-l-m, the other with wavenumber -m. This product equals 
the sum of two waves, one with wavenumber 2m and the 
other with wavenumber 0, which contributes to the mean 
flow. 

The mean zonal flow maintained by internal gravity 
waves is typically depth dependent, i.e., a shear flow. This 
"anti-diffusive" nature of gravity waves can be explained 
heuristically. Imagine a prograde (m>0) wave and a retro- 
grade wave (m<0) excited at the same radius. If the an- 
gular velocity, f2(r), of the medium increases inward, the 
prograde wave is Doppler shifted to lower frequency as it 
spirals inward {measured from the mean flow), whereas the 
retrograde wave is shifted to higher frequency. These fre- 
quency shifts are relative to the frequency at which they 
were generated, a gen, where the angular velocity is Qge„. 
This Doppler shifted frequency is 



cr[r) = aa 



+ m['[lg 



n{r)]. 



(13) 



Since the dissipation rate, and therefore damping distance, 
are strongly frequency dependent, the prograde and retro- 
grade waves dissipate at different depths. A prograde wave 
transports positive angular momentum, whereas a retro- 
grade wave transports negative angular momentum. There- 
fore, where a prograde wave is dissipated the mean zonal 



flow is accelerated and where a retrograde wave is dissipated 
the mean flow is decelerated. In this way two waves excited 
at a given radius with the same frequency and wavenumber 
but spiralling in opposite directions inward can modify the 
radial gradient of angular velocity, causing it to migrate up- 
ward toward the source of the waves and periodically reverse 



the direction of the radial shear. Plumb ( 1977 1 showed that 



two waves are unstable and will produce a shear, even in the 
absence of an initial shear. 



3.1 Quasi-Linear Approach 

To investigate wave-driven flows in detail the full set of fluid 
equations should be solved explicitly, so that the HARS ap- 
pearing in (12) can be calculated directly. However, this is 
computationally demanding as it requires good spatial reso- 
lution for both the waves (small spatial scales) and the mean 
flows (large scales) and good temporal resolution for the 
waves (short time scales) and the mean flow (long scales) rl 
To avoid these difficulties researchers employ a variety of ap- 
proximations, whereby the horizontally averaged Reynolds 
stress of a wave is evaluated according to linear theory. 

In the quasi-linear formalism the equation for the mean 
flow, such as that seen in equation (12), is subtracted from 
the full fluid equation to yield equations for the distur- 
bances such as u' and w'. The "mean fleld" approximation 
is then employed | |Herring]|1963[ ) in which the combination 
{u'w')r— < u'w' >r (where the subscript " r" denotes deriva- 
tion with respect to radius) is assumed small compared to 
other terms in the equation. This renders the equations lin- 
earized thus allowing the disturbances to be written as sim- 
ple sirmsoidal functions of space and time. The horizontal 
average of these fluctuations is then identified as a wave 
momentum. 

This treatment is similar to that of Bretherton (1966) 
in which the wave energy equation can be reduced to a con- 
servation equation for the wave action density: 



dA 

In 



+ V 



AVg = 0. 



(14) 



Here the action density A is equal to E/w*, where E is the 
wave energy and oj* is the frequency in the inertial framelj 
The product of the wave action density and the vertical 
group velocity, Vg yields the wave momentum flux, F. This 
conservation equation is derived in the absence of dissipative 
mechanisms. In order to use this conservation equation in de- 
scribing mean flow evolution, dissipation has to be included. 
This is done in an ad hoc manner, yielding an equation of 
the form: 



dA 
~dt 



+ V -F ^ -F/L 



(15) 



in which L is meant to represent a damping length. Expo- 
nential decay (in space) of the wave flux is then retrieved 
by neglecting dA/dt. This omission requires that the wave 



^ Clearly other combinations could work, but this is the simplest. 
^ By "wave" we mean a non-axisymmetric oscillating flow repre- 
sented by a single spectral wave number, m, not equal to zero; 
whereas a "mean flow" is an axisymmetric zonal flow with a spec- 
tral wave number of zero. 



■* For example, the QBO has a period of approximately 28 
months, but is thought to be driven by waves with periods less 
than approximately three days (Baldwin et al.|2001|l 
^ Note that there is a correspondence between the action flux 
and the Reynolds stress. 
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action vary slowly, similar to the assumption of WKB the- 
ory. This formalism results in the equation that is typically 
solved: 



& Bretherton (19671 studied the wave-critical level interac- 



9(7 

dt ^-^ r dr 









- + 



(16) 



where t represents a damping optical depth (Kumar et al 



1999h. A sum has been included to account for the flux from 



multiple waves, since the formalism described above applies 
to individual, non-interacting wave packets. Both of the ap- 
proaches described above neglect wave- wave interactions. 

In the treatment described above the fluid is "weakly 
turbulent", in the sense that the nonlinear terms involving 
wave self interactions and their influence on the mean flow 
are retained, but those involving wave-wave interactions are 
neglected. This approximation is valid when the divergence 
of the Reynolds stress term, {u'w'),. is small compared with 
the inertial term du' /dt, this limit is often described by re- 
quiring a small Froude number (w'/NL), where w' repre- 
sents a typical wave velocity, L represents a typical length 
scale and N represents the Brunt- Vaisala frequency. Other 
nonlinear effects, such as critical levels (where u — rnU/r 
approaches zero) and internal reflection (where a tends to 
N) are also not described in the above formulation (Booker 
fc Bretherto"nl|1967[ ). 



3.2 Attenuation Mechanisms 

Eliassen & Palm (1960) showed that in the absence of damp- 
ing, forcing or critical layers, no momentum could be trans- 
ferred from the waves to the flow. Stated another way, di- 
vergence of the Reynolds stress (12) is needed to maintain 
mean flows. The primary ways in which waves can be atten- 
uated are by radiative dissipation or by critical levels. In- 
ternal gravity waves are, in essence, thermal perturbations 
and therefore, are subject to thermal diffusion. As a wave 
propagates vertically, the wave's amplitude is decreased due 
to radiative diffusion. This dissipation leads to a nonzero 
divergence of the Reynolds stress, hence leading to angular 
momentum transfer and acceleration of the mean flow (12). 
The radiative dissipation of the wave depends on the wave 
properties, such as the frequency and the wavenumber. One 
can define a radiative damping distance (Kumar et al.|1999l 



d = 



Kk^m 



(17) 



where a is the wave frequency, k is the thermal diffusiv- 
ity, k is the wavenumber (m/r) and N is the Brunt- Vaisala 
frequency. This equation says that the higher the wave fre- 
quency and the smaller the wavenumber the less rapidly the 
wave dissipates. That is, the higher the frequency the less 
time there is to diffuse and the larger the wavelength the 
smaller the spatial gradients that drive diffusion. Radiative 
dissipation occurs continually as the wave propagates; there- 
fore the resulting mean flow acceleration is very gradual. 

Another way in which waves are attenuated is by crit- 
ical level absorption. A critical level is defined as that ra- 
dius where the mean zonal velocity, U, equals the horizontal 
phase speed of the wave, Cph = cr/k. In reality a wave packet 
has a range of horizontal phase speeds; therefore the criti- 
cal level is more appropriately referred to as a critical layer. 
The inviscid WKB solution is singular at this radius. |Booker| 



tion relaxing the WKB approximation and showed that, in 
the limit of large Richardson number {Ri = N^ /{dU/dz)'^), 
the wave energy density is strongly attenuated at a critical 
level with nearly complete momentum transfer to the mean 
fiow. More recent studies (IWinters & D'Asaro||1989l 11994 1 



indicate that some wave energy is transmitted through the 
critical level and some is refiected, but approximately one- 
third of the wave energy accelerates the mean fiow. Mean 
fiow acceleration at a critical layer is thus rapid and rela- 
tively local. 



4 RESULTS 

4.1 Two Wave driven shear oscillations 

4-1.1 Model Dependencies 

We simulate a mean flow oscillation maintained by two 
waves excited at our top boundary: a prograde and retro- 
grade wave each with the same longitudinal wavenumber 
and frequency. There is no background rotation in this case. 
Figure 1 shows the resulting mean angular velocity as a func- 
tion of time and radius for model M2 (see Table 1). Red and 
blue represent prograde and retrograde mean angular veloc- 
ity, respectively. This figure illustrates the salient features of 
a gravity wave driven oscillating shear fiow: (1) periodicity 
and persistence, (2) propagation of the mean fiow pattern 
toward the wave source in time and (3) a doubly peaked 
shear layer at all times. All of the oscillating models listed 
in Table 1 exhibit this fundamental behaviour. 

Figure 2 shows snapshots of the flow evolution during 
one cycle. Figures 2a and 2b show the full disk evolution of 
vorticity (niT^O), whereas (c) and (d) show the mean angular 
velocity as a function of radius. One can see both the mean 
fiow and the wave dynamics in this figure. In (c), the sec- 
ondary mean fiow (i.e., the deeper one) is prograde and in 
(a) the waves in this secondary layer are propagating in the 
retrograde (i.e, clockwise) direction. Note that since these 
waves are driven at the top boundary the energy, which is 
transported by the group velocity, spirals inward. Therefore, 
since the phase velocity of internal gravity waves is perpen- 
dicular to the group velocity, the phase velocity has a dom- 
inantly outward component in addition to its prograde (i.e, 
counter-clockwise) or retrograde (i.e., clockwise) tilt. This 
phase velocity is easily seen in movies of these simulations. In 
the snapshot illustrated in Figure 2a the phase is propagat- 
ing upward with a retrograde tilt in the secondary, whereas 
in the primary (the upper layer) the phase is propagating 
upward with a prograde tilt. The opposite configuration is 
occurring in Figure 2b, which is later in time. 

The existence of predominantly retrograde waves where 
the mean flow is prograde occurs because the prograde waves 
have been dissipated in this region; thus giving rise to the 
prograde mean fiow. This observation also means that some 
of the energy in the retrograde waves is able to propagate 
through the retrograde primary in Figure 2a (this is dis- 
cussed further in 4.1.3). 

The mean flow oscillation depends on the wave prop- 
erties (wavenumber, frequency and amplitude), as well as 
the properties of the medium in which the waves propa- 
gate. A series of models were run, varying these parame- 
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ters to investigate the various dependenciesPl Table 1 lists 
the various two-wave models and their parameterslj These 
dependencies will be discussed only briefly here since they 
have been reported previously ( Plumb|[l977 Wedi & Smo- 
larkiewicz|2006| . 

As is discussed in 4.1.2, we find that the mean flow oscil- 
lation period is determined largely (although not entirely) by 
radiative dissipation. Therefore, properties of the oscillation 
such as the period and depth can be qualitatively under- 
stood in terms of the radiative dissipation of waves; waves 
which are dissipated more rapidly produce shallower mean 
flow oscillations with shorter periods. This can be under- 
stood in terms of the gravity wave phase speed: waves with 
higher phase speed propagate further in a diffusion time. 
Therefore, gravity waves with larger frequency or smaller 
wavenumber produce mean flow oscillations that penetrate 
deeper and have longer periods. Likewise, larger thermal 
and viscous diffusivities lead to mean flow oscillations with 
shorter periods and shallower depths. Whereas the energy 
dissipation of the waves is predominantly thermal, the re- 
sulting angular momentum transfer is mediated through vis- 
cosity (12); therefore viscous diffusivity qualitatively affects 
the oscillation period and depth in a way similar to that of 
thermal diffusivity. 

Using the simpliflcd quasi-linear approach discussed in 



4-1.2 Mean Flow Evolution 

The major difference in the resulting behaviour in these sim- 



Section 3, Kim & MacGregor (20011 show that the nature 



of the mean flow oscillations produced (beyond period and 
depth) depends sensitively on the viscous dissipation. They 
find that as the viscous diffusivity is decreased, the mean 
fiows progress from steady shear (no oscillation) , to periodic 
oscillations, to chaotic reversals. 

The results presented here generally agree with those 
of Kim & MacGregor ("20011 with one important excep- 
tion. If the viscous diffusivity of model M2 is increased to 
5 X lO^^cm^s"^ (a factor of five), keeping other parame- 
ters fixed, the mean flow reverts to a steady shear flow. As 
the viscous diffusivity is decreased (M8-M10, displayed in 
Table 1) periodic reversals develop and the peak mean an- 
gular velocity increases. However so far we have not been 
able to produce chaotic solutions in these simple two-wave 
models. Our simulations with viscous diffusivity less than 
lO^^cm^s"^ are not well resolved (keeping other parameters 
fixed). The same behaviour can be recovered, in the oppo- 
site sense, by varying the amplitude of the forcing (Table 
1, Models M14-M16). As the forcing amplitude is increased 
fiows progress from steady to periodic oscillations and the 
peak mean angular velocity increases. Note, we have run 
cases with much smaller driving amplitudes (down by 100) 
and diffusivities (each down by 100) than those listed in Ta- 
ble 1. However, in those models a mean fiow does develop 
slowly, but after 200 years of simulated time they have not 
reversed and has reached only half the critical amplitude. 



^ Note again, that tlie amplitude is in spectral space, Vr ~ 
mij}/rp. 

^ The models presented here represent a small subset of the mod- 
els attempted. This subset represents those values of the various 
parameters that could be reasonably resolved in a timely manner. 



ulations compared to those of Kim & MacGregor (20011 is 



the presence of critical layers. As described in section 3, crit- 
ical layers cannot occur in a quasi-linear formulation such 



as that used in Kim & MacGregor (20011. In the following 



discussion of critical layers we continue to refer to the "pri- 
mary" layer as the mean fiow that is closest to the wave 
source (the top in this case) and the '"secondary" as the 
mean fiow below the primary; this is illustrated in Figure 
3b. We also refer to mean fiows as being "critical" when the 
peak angular velocity meets or exceeds the horizontal phase 
speed of the driven wave and "sub-critical" when it does 
not. 

Figures 3 (a-e) show snapshots of the mean angular ve- 
locity during one cycle of model M12. The dashed lines rep- 
resent the horizontal phase speed of the driven wave relative 
to the rotating frame of reference; angular velocities with 
amplitudes greater than this are "critical" . After a reversal 
of the mean fiow, the primary shear layer maintains its maxi- 
mum amplitude while the secondary shear layer slowly grows 
and moves upward toward the source. For example, when the 
primary layer is retrograde, as in Figures 3 (a-c), retrograde 
waves, observed relative to the mean retrograde fiow at the 
top boundary, have a lower frequency and prograde waves 
have a higher frequency (equation 13). As depth increases, 
the retrograde fiow decreases and becomes prograde within 
the secondary layer. This favors the transmission of the ret- 
rograde waves and the radiative dissipation of the prograde 
waves, which makes the amplitude of the prograde secondary 
mean fiow grow. As prograde waves continually dissipate the 
shear, and thus frequency shift, becomes stronger, causing 
the peak of the secondary flow to move toward the source 
in time. The growth of the secondary layer is relatively slow 
and takes approximately 85% of the total oscillation period. 
When the secondary layer is accelerated to an angular ve- 
locity approaching "critical" the prograde waves are rapidly 
dissipated, transferring the bulk of their angular momentum 
to the mean flow above the critical layer and bringing about 
the reversal of the primary. The reversal takes only about 
15% of the total oscillation period. The rapid acceleration 
associated with a critical layer begins before the secondary 
flow reaches a critical value. This happens for two reasons. 
First, as mentioned above, a wave packet has a range of hor- 
izontal phase speeds and therefore there is a critical layer. 
Second, there are regions in which the local (nonaxisymmet- 
ric) longitudinal velocity is critical, despite the mean flow 
being sub-critical. 

Rapid dissipation above the critical layer and the as- 
sociated angular momentum transfer is predicted based on 
the theories described in 3.2. In our simulations we see that 
this rapid dissipation is mediated by nonlinear wave-wave 
interactions which transport wave energy from the driven 
wavenumber to higher harmonics. This transport can be 
seen in Figure 4, which shows the wave energy as a func- 
tion of radius and wavenumber during growth of the sec- 
ondary shear layer (a) and during reversal of the primary 
(b), for model M9. During the bulk of the cycle most of 
the energy is in the driven mode (m=15), with minimal en- 
ergy in higher modes. However, once the angular velocity 
approaches the critical value energy is transferred from the 
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Table 1. Model parameters, m is the horizontal wavcnumber, / is the wave frequency {cr/(2n)) and A is the forced 
amplitude of the streamfunction, i/i, at the top boundary, k and f are the thermal and viscous diffusivities, D is the 
depth of oscillation (the depth is measured as that point were the mean angular velocity drops significantly below 
lO"'^ rad/s) and P is the period of the oscillation in days. AV and CP represent the maximum angular velocity 
of the mean flow and horizontal phase speed of the driven wave (a/rn = 2m/7r), respectively. The asterisk on the 
angular velocities for model MIO and M16 indicates that the pro- and retrograde peak angular velocities were not 
equal as in the rest of the models. The value listed is an average. 



Model 



i{fi Hz) A(10"gm/(cms)) K(10"cm2/s) u{l0^^cm^/s) B{10^cm) P{days) AY{lO-^rad/s) CP{lO-^rad/s) 





10 


10 






15 


10 






20 


10 






30 


10 






15 


8 






15 


15 






15 


20 






15 


10 






15 


10 







15 


10 




1 


15 


10 




2 


15 


10 




3 


15 


10 




4 


15 


10 


2 


5 


15 


10 


0.5 


6 


15 


10 


0.2 







8.8 


544.0 


6.30 


6.28 






4.4 


86.80 


4.28 


4.18 






2.9 


25.40 


3.18 


3.14 






2.0 


2.900 


0.65 


2.09 






3.1 


28.90 


3.33 


3.35 






10.2 


463.0 


6.44 


6.28 






17.5 


1620. 


8.53 


8.37 




0.2 


4.9 


144.7 


5.20 


4.18 




0.5 


4.9 


119.2 


4.70 


4.18 




5.0 


NA 


steady 


1.30* 


4.18 


0.1 




6.3 


189.8 


4.20 


4.18 


0.5 




5.4 


113.4 


4.30 


4.18 


3.0 




2.9 


50.90 


3.90 


4.18 


1 




4.39 


26.60 


4.85 


4.18 


1 




5.36 


405.0 


3.75 


4.18 


1 




NA 


steady 


0.14* 


4.18 



driven wave to higher harmonics where it can be more eas- 
ily dissipated (because of the larger spatial gradients), thus 
bringing about rapid reversal of the primary. Note, because 
of the quadratic nonlinearity, the energy is transferred only 
to horizontal modes that are multiples of the driving mode 
(m=15). 

The quasi-linear description of wave driven mean flow 
oscillations omits wave-wave interactions and relies on vis- 
cous diffusion of the mean flow for the reversal. That is, the 
right hand side (rhs) of (12) must meet or exceed the wave 
forcing (second term on the left hand side (Ihs) of (12)). 
Figure 5 shows that this is not the case in our simulations. 
Even during a reversal (Figure 5c), the viscous force acting 
on the mean flow is lower in amplitude than the divergence 
of the nonlinear Reynolds stress. This is not to say that 
viscous dissipation does not play a role. Viscosity, in addi- 
tion to thermal diffusion, dissipates the waves and it acts on 
the mean flow. However, the nonlinear wave-wave momen- 
tum transfer discussed above, followed by the dissipation of 
waves dominates the dynamics of the mean flow, not the 
dissipation of the mean flow. 

Therefore, the presence of critical layers makes these 
models qualitatively different than those described by the 
classical quasi-linear theory. This qualitative difference may 
produce quantitative differences in the depth and period of 
the oscillation. While we have attempted models with sig- 
nificantly reduced amplitudes and diffusivities, we have not 
been able to produce an oscillation in 5 x 10^ s of integra- 
tion time. However, a shear flow with an amplitude approxi- 
mately half the critical value has developed and is still grow- 
ing, but the growth rate is very slow and we have concluded 
that resolving an oscillation under such conditions is not 
numerically reasonable. 



4-1.3 The IGW Critical layer Interaction 

As discussed in 3.2, it has been demonstrated in several 
previous studies that when a gravity wave is incident on 
a critical layer, rapid dissipation and subsequent angular 
momentum deposition occurs, resulting in the acceleration 
of the mean flow. We have discussed this process above 
as the mechanism that brings about rapid reversal of the 
primary shear flow. However, Figure 3 (a-c) clearly shows 
that although the primary mean flow is indeed critical to 
the driven wave, the mean flow is not accelerated, despite 
the constant bombardment of waves. Furthermore, as shown 
in Figure 2, despite a retrograde primary layer, retrograde 
waves are evident (indeed dominant) below in the secondary 
layer. Clearly, waves are propagating past this apparent crit- 
ical layer. 

To investigate this somewhat surprising behaviour we 
spectrally decomposed the wave motion for model M2 in 
frequency space in order to distinguish between prograde 
and retrograde waves for model M2. This was done during 
different times in the cycle representing either prograde or 
retrograde primary mean flows. The results are shown in 
Figure 6. When there is a prograde primary mean flow (top 
two plots), prograde wave energy is transferred to higher 
wavenumbers as discussed in 4.1.2. However, the wave en- 
ergy is additionally Doppler shifted to higher frequencies. 
Because the shear is so strong, this Doppler shift is substan- 
tial, thus producing waves with higher phase speeds despite 
the increased wavenumbers. Waves with phase speeds higher 
than the phase speed of the driven wave (represented by the 
dashed line in Figure 6) will not recognize the mean flow as 
a critical level and therefore, propagate freely into the inte- 
rioirl On the other hand, retrograde waves incident on the 



* The Doppler shift of waves at critical layers was recognized by 
Pritts ( 1982 I and termed "self-acceleration". 
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prograde primary do not experience a critical level, so there 
is very little energy found at higher harmonics. The same 
behaviour (of the opposite sense) is observed when there is 
a retrograde primary (bottom two figures). The net effect of 
an IGW-critical layer interaction depends crucially on the 
sense (prograde or retrograde) of both the waves and mean 
flow. 

As mentioned in 4.1.2, when the mean flow is sub- 
critical there are regions that have non-axisymmetric lon- 
gitudinal velocities in excess of the horizontal phase speed 
and so are locally super-critical. Likewise, when the mean 
flow is super-critical there are regions in which the local 
longitudinal velocity is sub-critical. Given the complexity of 
natural systems, it is likely that such non-uniformity occurs 
and should be considered when discussing the propagation 
of waves near critical layers. 

As a wave approaches a critical level its vertical 
wavenumber is reduced substantially. Therefore, numerical 
simulations of critical layers require fine radial resolution. In 
order to test that the resolution employed was adequate to 
resolve the decrease in vertical wavenumber, we compared 
simulations with radial resolution 2.6 times finer and found 
that the same nonlinear transfer near a critical layer occurs. 
Similarly, the period and depth of the mean flow oscillation 
were not affected by increased radial resolution. Therefore, 
the flow behaviour in the vicinity of a critical layer should 
be adequately resolved. 

4.2 Forcing a Wave Spectrum 

The above results are observed when we force only one pro- 
grade wave and one retrograde wave, with the same am- 
plitude, wavelength and frequency relative to the rotating 
frame. Natural phenomena are clearly not as simple. A time 
dependent spectrum of waves is excited in the solar interior 
and in the Earth's atmosphere. To understand the genera- 
tion of mean flow oscillations in nature we must understand 
how multiple waves interact. 

We conducted several experiments varying the number 
and make up (frequency, wavenumber and amplitude) of the 
driven waves. These models are listed in Table 2. In mod- 
els MW1-MW3 we keep the amplitude of every wave fixed, 
but vary the number of waves forced. For example, model 
MWl forces wavenumbers 10 and 15, each with frequencies 
10jj,Hz and l^^Hz. Therefore four prograde and four retro- 
grade waves are forced, each with amplitudes the same as 
model M15. If one considers wave energy to be oc v^ then 
models MWl, MW2 and MW3 have approximately 1.4, 1.5 
and 3.5 times more energy input than M15 (considering the 
dependence of velocity on wavenumber). We find that as we 
increase the number of waves, it is increasingly difficult to 
produce an oscillating shear flow. MWl produces an oscil- 
lation, but neither MW2 or MW3 result in oscillatory solu- 
tions. In fact, the peak angular velocities (over the time run) 
for MW2 and MW3 are approximately 1000 and 100 times 
smaller, respectively, than MWl, despite having increased 
energy input. 

These smaller mean flows for MW2 and MW3 suggests 
that the forcing (second term on the left hand side of (12)) 
of the mean flow is lower, since the viscous diffusivity is held 
constant. To investigate this we plot the HARS, < u'w' >, 
for models MW1-MW3 and M15 as a function of radius in 



Figure 7a. There are a couple of interesting features to note. 
First, the HARS does not increase with increased wave en- 
ergy. In fact, MW2 and MW3 have a lower HARS compared 
with model MWl at the top of the domain. Second, the 
radial profile of the HARS varies substantially. The mag- 
nitude of the HARS for model M15 drops by nearly eight 
orders of magnitude over the radius range shown, whereas 
the multiple wave models drop by only two-three orders of 
magnitude. This could be partially due to the wavenum- 
bers and frequencies chosen; MW2 forces a lower frequency 
and wavenumber implying longer dissipation length by (17). 
However, the difference is not entirely due to this effect as 
MW3 forces both lower and higher frequencies (than M15) 
and yet still has a significantly less steep slope. The differ- 
ence in slope has a drastic effect on the forcing of the mean 
flow since that forcing depends on the radial gradient of the 
HARS. It is clear that the behaviour of the HARS for an 
ensemble of waves is signiflcantly different from that of the 
sum of individual, non-interacting waves. 

As discussed in Section 3, in order to approximate ve- 
locity fluctuations as simple sinusoidal functions of space 
and time the sum (u'w')r— < u'w' >r is assumed small in 
comparison to other terms in the equation. In Figure 7b, we 
show the ratio of the Reynolds stress to the HARS at a ran- 
domly chosen azimuth. We see that when only one prograde 
and one retrograde wave are forced, u'w' ~< u'w' > and the 
approximation employed in the quasi-linear approach could 
be valid. However, as the number of waves is increased u'w' 
becomes signiflcantly larger than < u'w' >. In these cases 
in order to neglect the sum of these two quantities {u'w')r 
must be small in comparison to du /dt. Stated another way. 



the Froude number for the wave (w'/Nz) (Andrews & Mcln- 



tyre 19761, where z is the vertical scale associated with a 



wave and the other quantities have their typical meanings, 
must be small. Using velocity and vertical wavelength scales 
typical in figure 7b and a Brunt- Vaisala frequency of lO^** 
one can estimate the Froude number to be approximately 
0.1 indicating that nonlinear wave- wave interactions are not 
negligible. The Froude number is however, just an estimate 
of the relative importance of the inertial term and the di- 
vergence of the Reynolds stress. In these simulations we can 
directly measure this ratio, which we show in the following 
sections. 



4-2.1 Effect of Amplitude 

Nonlinear wave- wave interactions are not only important for 
large amplitude waves. Models MC1-MC3 have been run 
with amplitudes and diffusivities each reduced by a factor 
of 100. To compensate for the lower diffusivities we increased 
the numerical resolution by a factor of two in the horizon- 
tal direction and by a factor of 3.75 in the radial direction. 
As for the models described above, when multiple waves 
are driven no oscillation develops in the 40 years simulated. 
This is not entirely surprising, the low amplitude and dif- 
fusivities imply significantly longer oscillation periods (see 
Section 4.1.1). We cannot realistically run simulations long 
enough to rule out oscillatory behaviour. However, as in 
models MW2 and MW3, the nonlinear wave-wave interac- 
tions are non-negligible. Figure 8a shows the radial gradi- 
ent of the HARS (dashed lines) and the radial gradient of 
the Reynolds stress (solid lines) for models MC2 and MC3. 
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Table 2. Models with multiple waves driven. Columns are denoted as in Table 1. 



Model 


m 


f{M Hz) 


A (lO^'^gm/s) 


k(10" 


cm^ 1 s) 


u(W^crr?ls) 


MWl 


10,15 


10,15 


0.5 


1 




1 


MW2 


5,10,15 


5,10,15 


0.5 


1 




1 


MW3 


5,10,15,20 


5,10,15,20 


0.5 


1 




1 


MCI 


15 


10 


0.01 


0.01 




0.01 


MC2 


10,15 


10,15 


0.01 


0.01 




0.01 


MC3 


5,10,15 


5,10,15 


0.01 


0.01 




0.01 


MC4 


5,10,15 


5,10,15 


0.022 


0.01 




0.01 


MSF5 


5,10,15 


5,10,15 


see text 


0.01 




0.01 


MSF6 


5,10,15 


5,10,15 


see text 


0.01 




0.01 


MSK5 


5,10,15 


5,10,15 


see text 


0.01 




0.01 


MSK6 


5,10,15 


5,10,15 


see text 


0.01 




0.01 


MFl 


5,10,15 


5,10,15 


see text 


0.01 




0.01 


MF2 


5,10,15 


5,10,15 


see text 


0.01 




0.01 



Figure 8b shows the radial velocity perturbation, w', as a 
function of time at a radius equivalent to 0.697?© for the 
same models. There are two main features to note. First, 
the amplitude of the HARS is two orders of magnitude lower 
than the Reynolds stress (Figure 8a). This just represents 
the fact that it is more difficult for low amplitude waves to 
force a mean flow, given the same integration time. There- 
fore, in order to neglect the nonlinear term, {u'w')r, one 
needs to show that this term is small compared to the in- 
ertial term, dw' /dt. Upon inspection of Figure 8b one can 
easily estimate the amplitude of the inertial term as approx- 
imately 10~^ — 10"'' depending on the model. The Reynolds 
stress at this radius is also approximately 10"'^; at best the 
Reynolds stress is 10% of the inertial term, hardly negligible. 
The amplitudes of these waves, depending on the wavenum- 
ber, varies between 50 cm/s and 150 cm/s. These amplitudes 
combined with those modeled in MW1-MW3 represent any 
reasonable estimates for the wave amplitudes in the solar 
tachocline and radiative interior. 



4.2.2 Effect of Spectra 

The models listed above force waves with constant stream- 
function, i/), which makes the radial velocities proportional 
to m. This gives a wave energy "flux" per mode proportional 
to m'^. We ran five models in addition to models MC1-MC3 
to investigate the dependence of the results on the wave 
spectrum. These models are listed in Table 2 as models 
MC4 representing a model with wave flux proportional to 
m^ but with an energy input larger than that of MC3 by 
a factor of five. Models MSF5 and MSF6 represent models 
with a wave fiux which is fiat in real space (constant wave 
energy fiux per mode), whereas models MSK5 and MSK6 
represent models with a wave fiux proportional to m^/~^, 
similar to that proposed by Kumar et al. (19991. Since the 



amplitudes of the individual waves in these models varies 
substantially, we do not list them in Table 2. Furthermore, 
models appended with a "6" represent models with a total 
wave energy input five times larger than models appended 
with a " 5" . The amplitudes of the velocities for models ap- 
pended with "5" vary from 20 cm/s to 150 cm/s, whereas 
those appended with "6" vary from 60 cm/s to 340 cm/s. 
These models were run for a total of 200 years. Although, 
none showed oscillatory behaviour in that time, the models 



with higher input energy (MC4, MSF6, MSK6) do maintain 
strong shear fiows. Our results do not rule out oscillations 
with periods longer than hundreds of years. 

The results for these models are shown in Figure 
9 which represents the same physical quantities as por- 
trayed in Figure 8: (a) and (b) represent the lower in- 
put energy models (MC3,MSF5,MSK5), while (c) and (d) 
show the various spectra models with higher input energy 
(MC4,MSF6,MSK6). There are several interesting features. 
First, the amplitude of the HARS, represented by dotted 
lines, is approximately four orders of magnitude smaller than 
the non-averaged Reynolds stress for the lower energy mod- 
els. Therefore, {u'w'),.— < u'w' >,.~ {u'w')r and in order 
to neglect these terms {u'w')r must be small compared to 
dw' /dt. Upon inspection of Figure 9b one sees the amplitude 
of the inertial term is approximately 4 x 10""^ (this ampli- 
tude varies very little with spectrum); examining Figure 9a 
at the same radius (0.697?©), one sees that the divergence 
of the Reynolds stress, {u'w')r, has virtually the same am- 
plitude. Therefore, the Reynolds stress cannot be neglected 
and the problem is fundamentally nonlinear. Finally, one can 
see that the slopes for models MSF5 and MSK5 (represented 
by blue and red, respectively) are significantly different than 
the slope for model MC3 (represented by black), despite 
having the same frequencies and wavenumbers driven. The 
radiative wave damping for an individual wave, described by 
(17), implies that the slope should depend only on frequency, 
wavenumber, diffusivity and Brunt- Vaisala frequency. These 
parameters do not change between models MC3, MSK5 and 
MSF5, yet the slopes change substantially. This observation 
implies that not only are wave- wave interactions relevant but 
that they are spectrum dependent. The higher input energy 
models MC4, MSK6 and MSF6 show virtually the same be- 
haviour, with one exception. In these models mean fiows are 
more readily generated over the same period. Whereas the 
HARS is significantly larger than in the lower energy mod- 
els, it is still lower than the un-averaged Reynolds stress by 
one to two orders of magnitude. One can make the same 
comparison between {u'w')r and dw' /dt and find the same 
result; namely that the nonlinear terms are non- negligible. 

It is also worth pointing out that the wave fiux predicted 



by Kumar et al. ( 1999 1, represented by the red lines in Fig- 



ure 9 appear to be more efficient at driving mean fiows. In 
both the low energy models and the high energy models the 
HARS is closer to the non-averaged Reynolds stress than in 
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the other two models. The wave-wave interactions appear 
to be spectrum dependent. However, in addition, it appears 
that the (in)coherence of these interactions and their ability 
to force a mean flow is spectrum dependent. 



4-2.3 Effect of Initial Conditions 

All of the models described above start with no mean flow. 
Therefore, if the waves cannot force a mean flow, by con- 
struct the problem is non-linear, in the sense that one can- 
not separate the flow into a mean (with large amplitude) 
and a wave (with smaller amplitude). Therefore it may not 
be surprising that these models exhibit such nonlinear be- 
haviour. To investigate whether wave-wave interactions are 
influenced by the presence of a mean flow we ran two addi- 
tional models, labeled in Table 2 as MFl and MF2. These 
models have a mean flow forced from the outset: MFl has 
a sinusoidal function forced in the top 10^ cm (or 2% of the 
radius of the solar interior), while MF2 has a half-sinusoid 
function in the top 10® cm. In both models the mean flow 
is specified to have a maximum amplitude of IC'cm/s and 
is forced with a spectrum of waves which all contribute the 
same amplitude to v^, making the total Vr peak at 100 cm/s. 
The mean flow is enforced such that wave-wave interactions 
and wave-mean flow interactions are allowed through the 
nonlinear terms. The results of these two models are shown 
in Figure 10; similar to figures 8 and 9 the upper panel 
shows the radial derivative of the Reynolds stress (solid 
lines), and the radial derivative of the HARS (dotted lines). 
The Reynolds stress depicted is due to only wave-wave in- 
teractions, that is, wave interactions with the imposed mean 
flow is not included in the value shown in Figure 10. Black 
represents model MFl, while blue represents MF2. In these 
models the HARS is still an order of magnitude lower than 
the non-averaged Reynolds stress (Figure 10a), but again 
the Reynolds stress has an amplitude similar to dw' /dt as 
measured in Figure 10b. Even in models in which a mean 
flow is prescribed and the wave amplitude is prescribed to 
be 1% of the mean flow amplitude, nonlinear interactions 
are non-negligible. 



5 DISCUSSION 

We conducted highly simplified, but fully nonlinear, numer- 
ical experiments of internal gravity waves in the solar ra- 
diative interior to understand whether such waves can force 
oscillating shear flows. Our results are threefold. First, we 
produce oscillating shear flows for models in which a single 
horizontal wavelength is forced as equal amplitude prograde 
and retrograde propagating waves. However, most of these 
highly simplified models exhibit critical layers which domi- 
nate the mean fiow dynamics. Second, when waves are inci- 
dent on a critical layer the amount of energy deposited com- 
pared to the amount of energy transmitted varies depending 
on the sense of the mean fiow those waves propagate into. 

Finally, when a spectrum of waves is forced the standard 
quasi linear approximation which neglects wave-wave inter- 
actions is inadequate. In fact, none of the models we con- 
ducted could be adequately described by quasi-linear theory, 
either because of the presence of critical layers or because of 
the dominance of wave-wave interactions. Such wave-wave 



interactions are non-negligible from waves ranging in am- 
plitude from 50-1000 cm/s and diffusivities ranging from 
10® — lO^^cm^ /s for various spectra of driven waves, with or 
without an enforced mean flow. This range covers any rea- 
sonable estimate for flow velocities and turbulent diffusivi- 
ties in the solar tachocline. Lower diffusivities would make 
the nonlinear wave-wave interactions even more important. 
Therefore, it appears that in order for quasi-linear theory 
to be valid for flows in the solar tachocline one would need 
unreasonably low amplitude waves combined with unreason- 
ably high diffusivities. 

The ability for waves to enforce uniform rotation in the 
radiative interior depends sensitively on the ability to set 
up an oscillating shear layer in the solar tachocline. It has 
been shown previously that that ability depends sensitively 
on the spectrum and amplitudes of the waves driven. Here, 
we further show that the simplified theory used to describe 
those flows is inadequate for various spectra, amplitudes and 
diffusivities. 
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Figure 1. Mean angular velocity as a function of radius and time 
for model M2. Only the outer 30% in radius of our simulated sta- 
ble region is depicted here. Prograde motion is represented by red, 
while retrograde motion is represented by blue. The main features 
of a gravity wave driven mean flow oscillation are recovered. 
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Figure 2. Flow dynamics over a cycle for model M2. Vorticity 
is depicted in (a) and at a later time in (b). Positive vorticity is 
represented by yellow/white and negative vorticity represented by 
dark (blue/black) colors, (c) and (d) represent the mean angular 
velocity at the times corresponding to (a) and (b), respectively. 
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Figure 3. Mean angular velocity at several times over one cy- 
cle for model M9. Dashed lines represent the phase speed of the 
driven wave. 
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Figure 4. Wave energy as a function of wavenumber and radius 
for model M9. During growth of the secondary flow (left) wave 
energy is found primarily in the driven wavenumber (m=15). Dur- 
ing reversal of the primary flow (right) wave energy is found dis- 
tributed amongst all possible modes. The transfer of energy to 
higher harmonics allows rapid wave dissipation and associated 
angular momentum transfer. 
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Figure 5. One half cycle of model Mil. The solid line is the 
divergence of the nonlinear Reynolds stress (the second term LHS 
of equation 12, and the dotted line is the divergence of the viscous 
stress acting on the mean flow (the RHS of equation 12), both 
in cm/s^. The nonlinear term is always larger than the viscous 
term acting on the mean flow. 
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Figure 6. Wavenumber (m)-Frequency (f) diagram of wave en- 
ergy for prograde and retrograde waves at O.69R0 for model M2. 
Top panels show this when there is a prograde primary mean flow 
at this radius; whereas the bottom panels are when a retrograde 
primary mean flow exists. The diagonal line represents the mean 
angular velocity. Wave energy below this line should be trapped 
by a critical layer. This shows that waves of the same sense as 
the mean flow are doppler shifted to higher frequencies, so that 
significant wave energy is able to propagate past the critical layer. 
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Figure 7. (a) Horizontally averaged Reynolds stress (HARS) ver- 
sus radius for models with multiple waves forced, compared with 
models with only a single (prograde and retrograde) wave forced 
(M15). (b) Ratio of Reynolds stress at a randomly picked azimuth 
to the HARS. The linetype is the same as in (a). As multiple waves 
are forced, the ratio strays increasingly from one. 
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Figure 8. (a) Radial derivative of the HARS (dotted lines) com- 
pared with the radial derivative of the non-averaged Reynolds 
stress (solid lines) with black representing model MC2 and blue 
representing model MC3, (b) radial velocity perturbation as a 
function of time at 0.69Rq and a randomly chosen azimuth. In 
both cases the inertial term is the same order as the nonlinear 
terms 
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Figure 9. Similar to figure 8, (a) shows radial derivative of the 
HARS (dotted lines) compared to the non-averaged Reynolds 
stress term (solid lines) for models MC3 (black), MSK5 (red) 
and MSF5 (blue); (b) radial velocity fluctuation as a function of 
time at 0.69 Rq and a randomly chosen azimuth for the models 
portrayed in (a); (c) similar to (a) but for models MC4 (black), 
MSK6 (red) and MSF6 (blue); these models have a total energy 
input five times larger than those shown in (a); (d) same as (b) 
but for the models portrayed in (c). In all cases the inertial term 
is the same order as the nonlinear terms. 
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Figure 10. Similar to Figures 8 and 9. (a) Radial derivative of 
the HARS (dotted lines) compared to the non-averaged Reynolds 
stress term (solid lines) for models MFl (black) and MF2 (blue), 
(b) Radial velocity perturbations as a function of time at 0.69-Rq 
and a randomly chosen azimuth. In both cases the non-averaged 
Reynolds stress is larger than the HARS and approximately equal 
to the inertial term. 
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